Higher-order Time Integration of Stochastic Differential Equations and Application to Coulomb Collisions

We have examined the extension of Langevin-equation Monte-Carlo algorithms for Coulomb collisions from the conventional Euler O(∆t)-order time integration to the next higher [Milstein-O(∆t)] order of accuracy. This examination has yielded many results, some of interest specifically for the Coulomb-collision problem, some of fundamental importance for integration of stochastic motion on a spherical (and, more generally, any curved) surface, and some of considerable utility for integration of multi-dimensional systems of stochastic differential equations (SDE’s). In one common Langevin-equation approach, the angular scattering is treated with a combination of near-Cartesian stochastic velocity-direction kicks, in a unit-vector frame that is rotated so that at the beginning of each timestep the third axis is aligned with the velocity direction. We find that in such schemes, the angular component of the collisional scattering cannot be extended beyond the Euler order. Instead, the extension to higher order proceeds via a formulation of the angular scattering directly as SDE’s in the two fixed-frame spherical-coordinate variables. This extension has been implemented, and results have been obtained for Coulomb collisions showing the improved [O(∆t) vs. O(∆t)] convergence. Our algorithm requires generation of stochastic “area integrals” that represent higher-order coupling between these variables, in addition to generation of the Wiener-process increments. We have found a simple but accurate approximation to the joint probability density function (PDF) of the Levy-area and the Gaussian Wiener-process increments that considerably simplifies the sampling of the area integrals compared to the direct method of Gaines and Lyons [SIAM J. Appl. Math. 54, 1132-1146 (1994)], reduces the dimensionality of the sampling function from 2 to 1, and allows for a direct implementation using simple analytical formulas as an alternative to tables. We have also developed a direct approach to the conditional sampling problem involved in adaptive integration of SDE’s. This may be computationally more efficient than the quadrature-formula approach of Gaines and Lyons [SIAM J. Appl. Math. 57, 1455-1484 (1997)], because it avoids generation and storage of, and computations involving the many random numbers that representing the Wienerprocess increments over time intervals much finer than the actual numerical time step at a given level of refinement. Our approximation to the Levy-area-Wiener-displacement PDF is very helpful here because its use reduces the dimensionality of the conditional sampling function from 4 to 3. The computational memory required for storage of the associated table is then reduced from large to quite modest, and for many SDE calculations will not represent a large component of their memory load. In our direct explicit Milstein implementation for the Coulomb-collision problem (and the associated simpler spherical diffusion problem) the equilibrium state is determined by a dynamic balance of drag and scattering terms. Small numerical errors in either may cause significant errors in the computed equilibrium. Two-step algorithms are therefore considered, in which a predictor step is used to estimate centered (Stratonovich-calculus) or forward values of the stochastic variables. These may be more accurate because they can eliminate the need for the dynamic balance in the determination of the equilibrium state. * Work performed for US DOE by LLNL under Contract DE-AC52-07NA27344 and by UCLA under grant DE-FG02-05ER-25710.